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Abstract. We present an efficient technique to study the ID evolution of instability-generated structure in winds 
of hot stars out to very large distances (~ 1000 stellar radii). This technique makes use of our previous finding 
that external forces play little role in the outer evolution of structure. Rather than evolving the entire wind, as 
is traditionally done, the technique focuses on a representative portion of the structure and follows it as it moves 
out with the flow. This requires the problem to be formulated in a moving reference frame. The lack of Galilean 
invariance of the spherical equations of hydrodynamics is circumvented by recasting them in a pseudo-planar form. 
By applying the technique to a number of problems we show that it is fast and accurate, and has considerable 
conceptual advantages. It is particularly useful to test the dependence of solutions on the Galilean frame in which 
they were obtained. As an illustration, we show that, in a one- dimensional approximation, the wind can remain 
structured out to distances of more than 1300 stellar radii from the central star. 
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1. Introduction 

The line-driven stellar winds of hot, luminous, OB-type 
stars are subject to a strong line-deshadowing insta- 
bility (e.g. Owocki & Rybicki HUgl Feldmeier |2TjUT|) . 
Hydrodynamical studies aimed at following the nonlinear 
evolution of this deshadowing instability (Owocki, Castor 
& RybickiUHEEl Feldmeier et al lTW71 Runacres & Owocki 
2002 : - hereafter Paper I) show that within a few stellar 
radii of the surface the flow becomes highly structured, 
with gas concentrated in dense clumps, and pervaded by 
strong shocks. To determine the full physical and observa- 
tional significance of such structure, it is important to un- 
derstand its subsequent development at scales beyond the 
few stellar radii covered in such radiation- hydrodynamical 
simulations of its initial formation. For example, for a star 
such as £ Pup, nearly half of the observed thermal radio 
flux is understood to originate at distances of more than 
100 i?*. For non-thermal emitters such as Cyg OB2 No. 9, 
shocks are needed beyond 500 i?* (Van Loo et al. I2UU4|I . 
For some stars (such as £ Pup) it has been suggested 
that a significant contribution to the X-ray flux originates 
beyond 100 i?* (Hillier et al 1993). Recently, however, 
Kramer et al. i|20Q."i) have found fits to X-ray emission 
lines that indicate a formation region much closer to the 
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star (r ^ 5 i?»). Finally, Grosdidier et al. I|1998fl have sug- 
gested that some of the structure found in ring-nebulae 
around Wolf-Rayet stars might be an imprint of stellar 
wind structure. 

These considerations demonstrate the importance of 
modelling the dynamical evolution of instability-generated 
structure out to very large distance scales of order ~ 
1000 i?*. Because of the computational expense of the 
non-local integrals for calculating the radiative force cen- 
tral to the developing instability, full radiation hydrody- 
namical models have generally been limited to the dis- 
tances of order ~ 10 i?*. But in Paper I we showed that 
such radiative forces become of negligible importance be- 
yond distance of a few times ~ 10 i?*, and so within a 
hybrid approach that switches to a much less costly, pure- 
hydrodynamical model of such outer regions, we were able 
to extend high-resolution, ID instability models to dis- 
tances of order ~ 100 Unfortunately, even with this 
approach, extension to still larger scales again quickly be- 
comes quite computationally expensive, since it requires 
both a large number of depth points and a evolution of the 
model for the extended time required for advected struc- 
ture to relax over such scales. 

To address this problem, the present paper introduces 
a new pseudo-planar, moving periodic-box approach. This 
further reduces the complexity of the models by trans- 
forming the analysis into the mean rest frame of some 
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representative section - the box - of the expanding flow 
structure. Then, to account approximately for the sec- 
ular radial changes associated with spherical expansion, 
the flow variables and their governing equations are re- 
cast in forms that resemble those for simple planar flow, 
the "pseudo-planar" form. Finally, under the assumption 
that the chosen section represents a randomly character- 
istic sample of the flow structure, its internal evolution is 
then isolated by assuming periodic boundary conditions 
to link the inner and outer edges of the box. 

In the following we first review the radiatively 
driven models that serve as basic input to our approach, 
and then (fJ2| formally introduce and develop the pseudo- 
planar, periodic box formalism. We next (<2J apply this 
to a number of test problems, and ffl5)l compare the re- 
sults with those of traditional, radiatively driven models. 
In both sections, we devote particular attention to the 
effect of advection and show that the ability to test the 
dependence of problems on the Galilean frame is a major 
advantage of the periodic box model. After illustrating 
(§Hj the application of the method in a model with struc- 
ture extending to 1300 i?*, we conclude (J7|) by summa- 
rizing the advantages, potential applications, and future 
extensions of a periodic-box approach. 

2. Radiatively driven models of outer wind 
structure 

In Paper I we investigated structure up to a distance of 
100 i?*, using Eulerian, one-dimensional, time-dependent 
hydrodynamical models that take into account the insta- 
bility of the driving. In these models, the material is com- 
pressed into a sequence of narrow, dense shells, bounded 
by shocks. These shells expand at roughly the sound speed 
as they move out at approximately the terminal velocity 
of the wind. There are supersonic velocity differences be- 
tween individual shells, causing them to collide and form 
new shells. The importance of similar shell-shell collisions 
for the production of X-rays has already been pointed out 
by Feldmeier et al. I|1997[) . We found that these collisions 
effectively hinder the decay of the structure initiated in the 
inner wind, so that the clumps can survive to substantial 
100 i?*) distances. 

In modelling the distant wind structure, it is necessary 
to maintain a relatively fine grid spacing to resolve the 
often quite narrow dense clumps. For a radiatively driven 
calculation, we use a grid of 10 637 points. The initial 1025 
points have a spacing that increases linearly from 0.001 to 
0.01 i?* over the range 1 — 5 i?*. Beyond 5 i?* we use a 
constant grid spacing of 0.01 i?*. This is also the spacing 
of the box models. 

Another important factor in the outer evolution of 
structure is the energy balance in the wind. The gas in 
our simulations cools both adiabatically and radiatively. 
The combined effect of this cooling is balanced by pho- 
toionisation heating from the star's ultraviolet radiation. 
We mimic the effect of radiative heating by setting a floor 
temperature, below which the temperature is not allowed 



to drop. The value of this floor temperature influences 
the expansion speed of the shells. For a low value of the 
floor temperature, the dissipation of the structure will be 
slower. In Paper I, the floor temperature was taken equal 
to the effective temperature. This is a crude approxima- 
tion and probably results in a wind that is too warm. As 
a first step towards a more realistic treatment of the en- 
ergy balance, we used the detailed ionisation and thermal 
equilibrium models by Drew 1)1989(1 . These take into ac- 
count the cooling by heavy element lines and predict an 
outward decreasing temperature that quickly reaches val- 
ues substantially below the effective temperature. These 
temperature profiles can be adequately approximated by 
the expression 



(Bunn & Drew rH)92[) . which we used in the current mod- 
els, with the velocity assumed to follow the usual "beta- 
law" form v(r) — Woo(l — R^/r) 13 with exponent (3 — 0.7. 
Note that the temperature profile levels off at 0.28 T c g. 
Although this represents a modest improvement over sim- 
ply taking the effective temperature, it is still quite unre- 
alistic. The Drew models were only calculated to 10 
Most importantly, they assume a smooth outflow. There 
are numerous ways in which the inhomogcncity of the out- 
flow could influence the energy balance. A detailed inves- 
tigation of the energy balance in a structured wind, how- 
ever, is beyond the present scope. 

Finally, the amount of structure initiated in the inner 
wind strongly depends on the line driving parameters, in 
particular on the cut-off parameter K max that limits the 
maximum line strength (Owocki, Castor & Rybicki 1988). 
For purely computational reasons, this parameter is usu- 
ally set to artificially low values. We found that, with the 
relatively fine resolution of our calculations, it is possible 
to set this parameter to less artificial values. In the current 
paper, we have used K max — O.Iko, where the opacity con- 
stant kq is related to the actual strength of the strongest 
line. This is a factor of hundred larger than in Paper I. The 
effect is to include a number of strong lines that become 
optically thin only for very large velocity gradients. This 
allows for much stronger rarefactions and shocks. The re- 
sulting models are extremely structured, as can be seen 
from the snapshot (Fig. f^) . 

In summary, we can say that the amount of clump- 
ing in our simulations is largely determined by the grid 
spacing, the floor temperature and K max . On the other 
hand, we find that the clumping does not depend on the 
radiative force beyond ~ 30 i?* (Paper I). This reduces 
the outer evolution of structure to a pure gasdynamical 
model. As the evaluation of the radiative force dominates 
the computation time, this allows us to construct vastly 
more economical models, which we present in the following 
section. 
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Fig. 1. Velocity, density and temperature for a radiatively 
driven model. The dashed lines in the two upper panels 
represent the time-averaged variable. 

3. A pseudo-planar moving periodic box formalism 

3.1. Motivation and outline 

Even without the evaluation of the radiative force, 
the modelling of structure out to very large distances 
(~ 1000 i?*) is still expensive and the fine grid spacing 
needed to resolve the shells results in impractically large 
files. While such models are marginally feasible in ID, 
generalizing them to higher dimensions would be compu- 
tationally prohibitive. In this section, we present a tech- 
nique that does not suffer from this disadvantage, and has 
considerable conceptual benefits. 

First, let it be noted that the structure generated by 
the instability, apart from being stochastic, is also quasi- 
regular in the sense that similar features are repeated over 
time. Therefore it is not necessary to keep track of the 
whole stellar wind during the duration of the simulation. 
It is enough to select a limited but representative portion 
of the structure, and follow this as it moves out with the 
flow. Following a portion of the wind entails transforming 
the conservation equations to a moving reference frame. 
This is not possible directly, as the spherical equations of 
hydrodynamics are not invariant under a Galilean trans- 
formation. This problem can be circumvented by rewriting 
the equations in a pseudo-planar form (see below). 

A fundamental point in the present analysis is that the 
timc-dcpcndcnt hydrodynamical variables can be viewed 
as having two distinct components, varying on different 



scales. The first is a rapidly oscillating variable, varying on 
short spatial and temporal scales. The second is a slowly 
varying function of radius, corresponding to the secular 
evolution of the time-averaged variable. As an example, 
consider the density in Fig. The gas in this simulation is 
highly clumped, with density variations of several orders 
of magnitude over less than a stellar radius. The time- 
averaged density however, is a much more well-behaved 
function, decreasing steadily outward. Indeed, one of the 
more reassuring properties of the models including the de- 
shadowing instability is that the mean variables closely re- 
semble the results from CAK theory. In our pseudo-planar 
reformulation of the hydrodynamical equations we try to 
absorb, as much as possible, the secular evolution in the 
scaling of the variables. 

3.2. Pseudo-planar scaling 

For the sake of clarity, let us recapitulate the spherical 
equations of hydrodynamics. In Eulerian form, the one- 
dimensional (ID) time-dependent equations for conserva- 
tion of mass, momentum, and energy are: 



dp 

at 

dt 

de 

at 



1 d(r 2 pv) 
r 2 dr 
1 d(r 2 pv 2 ) 
r 2 dr 
d(r 2 ev) 
dr 



= 



1 



-5 + / 
dr 

p d{r 2 v) 
r 2 dr 



Q. 



(2) 
(3) 
(4) 



Here / is the sum of the external forces (gravity and radia- 
tive driving) acting on the gas and Q the power emitted 
by radiation per unit volume. The internal energy density 
e (in erg/cm 3 ) is related to the pressure by the perfect 
gas law p = (7 — l)e = pkT / (/ima) , which supplements 
the set of equations. The other symbols have their usual 
meaning. In this paper, we assume a perfect monatomic 
gas and use 7 = 5/3 , unless otherwise specified. 

The first step in our reformulation of the equations of 
hydrodynamics is to scale out, as much as possible, the 
secular components of the density, velocity and pressure. 
From the continuity equation it can be seen that the sec- 
ular expansion of the gas at large distances from the star 
(where we have reached the terminal velocity) causes the 
density to fall of as l/r 2 . This suggests the following def- 
inition of the scaled density p: 



2 - 2 
r p = r p, 



(5) 



where ro is a fiducial radius. As the mean temperature 
in the present models is almost constant, it is logical to 
scale the pressure in the same manner as the density. The 
scaled internal energy then follows from the perfect gas 
law p = (7 — l)e. The mean velocity in the outer wind is 
roughly constant and need not be scaled. 

Rewriting Eqs. m terms of the scaled variables, 

we obtain the following conservation equations: 



dp 

dt 



d (pv) 
dr 







(6) 
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dt dr 

de d (ev 
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dt dr 



or r 



_dv 
^ dr 



Q 



2pv 



(7) 
(8) 



where the external force and the heating rate have been 
scaled in the same way as the density and the pressure. 
These equations are pseudo-planar in the sense that they 
describe the spherical problem but formally resemble the 
planar equations of hydrodynamics. The only formal dif- 
ference with the planar equations are the geometric source 
terms appearing in the momentum and energy equations. 
These express the secular evolution of momentum and en- 
ergy in spherical geometry. The outer wind is essentially 
expanding at a constant velocity while the temperature is 
kept constant by the balance between heating and cooling. 
The geometric source term 2p/r in the momentum equa- 
tion represents the pressure gradient associated with this 
secular expansion. The geometric source term 2pv/r in the 
energy equation expresses the secular adiabatic cooling of 
the gas. 

3.3. Transformation to moving box 

Let us now transform the pseudo-planar equations from 
the stellar rest frame to a frame that is moving at a con- 
stant velocity v . We refer to this moving frame as the 
box. We define the position coordinate x within the box 
by the Galilean transformation 



r = R + v a t + x, 



(9) 



where Rq is the initial inner radius of the box. For any vari- 
able a depending on position and time we can then convert 
from a(r,t) to a(x,t). We denote the time-derivative at 
constant x by d/dt. We then have 



da 
"dT 



da 
~dt 



da 
~dt 



+ Oo • V) a 



(10) 



If we further introduce the velocity w with respect to the 
box velocity, i.e. v = vq + w, then we can rewrite the 
equations of hydrodynamics as 



dp 
dt 
d (pw) 



d (pw) 

dx 
(pw 2 ) 







d 



dt 



dx 
d (ew) 



de 

dt dx 



dp + , | 2p 

dx Rq + vot + x 

„dw ~ 2p(v + w) 

~p~a~ + Q- -5—, — r~; — 

dx Rq + Vot + x 



(11) 
(12) 

(13) 



Finally, we impose periodic boundary conditions on 
the box. The box method has a number of advantages that 
should be clear already. Not only is the number of grid 
points reduced with respect to a traditional calculation 
where the whole wind is evolved, it is also not necessary 
to increase the number of depth points if the structure 
needs to be followed further. The computing time is less, 
due to the smaller number of grid points. Furthermore, 
any explicit time-stepping method is subject to a stability 



condition on the Courant number C — cAt/Ax, where c 
is the maximum propagation speed of information, At the 
time step and Ax the grid spacing. In order to be stable, C 
must be smaller than one (this is the Courant-Friedrichs- 
Lewy condition), and preferably a lot smaller. Due to the 
slower speed at which features are advected over the nu- 
merical grid, the condition on the Courant number is less 
restrictive for the box model. The gain is substantial: per- 
mitted time steps increase from ~ Ax/v^ to a few times 
Ax/a, where a is the sound speed. In the simulations pre- 
sented here, Voo/a ~ 150. 

3.4. Effect of periodic boundary conditions 

General principles As a consequence of the periodic bound- 
ary conditions, structure is always dragged along with the 
box. To ensure that the time scale over which the structure 
evolves is comparable to the time-scale associated with its 
outward movement, the speed of the box must be roughly 
the same as the terminal velocity. This is clear from the 
fact that it is obviously not possible to evolve structure 
cheaply by taking a very large box speed. Furthermore, 
as shell collisions are very important for the evolution of 
structure (Paper I), the box should be big enough to con- 
tain more than one shell, even at late times. Finally, the 
periodic boundary conditions move features such as dense 
shells from one side of the box to the other. The box 
method is inherently incapable of describing variations on 
scales beyond a box length. This means that the length of 
the box should be small compared to the distance covered 
by the box. 

Restriction of box size from energy equation A quantitative 
restriction on the box size can be derived from the secular 
change of the energy density over the length of the box. 
Using the energy equation (|13f) we can write the relative 
change in the secular component of the energy density 
over a box length L as 



1 



2p(v + w) 
Ro + vot + x 



dt. 



(14) 



where r = L/vq. If we assume that the change in p is slow 
compared to 1/t this gives 



2(7- 1) In 1 + 



L 



Ro 



(15) 



Using x — L/2 and Ro = 50i?» the requirement e c < 0.25 
gives L < 11.5.R*. 

Restriction from clumping factor A different measure of the 
effect of periodic boundary conditions is given by their 
influence on the clumping factor. When a shell crosses the 
boundary of the box, it effectively jumps a distance L in 
the stellar rest frame. If the clumping factor changes over 
length scales ~ L, this introduces an error, which we can 
estimate in the following simplified model. Let us assume 
a single shell in the box, with an expansion speed a c . If 
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the shell has a width Zl at the left side of the box then 
it has a width Zl + Al at the right side of the box, where 
Al = deL/voo. Note that the relevant time-scale for the 
expansion is not the time needed for the shell to cross the 
box, but the much shorter time L/vq needed for the shell 
to move a distance L in the stellar rest frame. We can then 
write the relative change in the clumping factor caused by 
applying the periodic boundary conditions as 



£cl 



/cl,L — /cl.R 



cl,R 



0; 



L 



(16) 



where we have used the fact that the clumping factor can 
be approximated by the inverse of the volume filling fac- 
tor, which in turn is the fractional width occupied by the 
shell (see e.g. Paper I). This equation shows that the error 
caused by applying periodic boundary conditions is large 
if the shells expand fast, the outflow velocity is small, the 
box is big or the shells are narrow. This last possibility is 
explained by the fact that for a given L,a e and Voo, the 
expansion Al is constant and therefore relatively more im- 
portant for narrow shells. From the models presented later 
in this paper (Sect.EJ), we can derive typical values to use 
in Eq. I|16|l . At w 170i?*, where the first shell crosses the 
box (Fig. 111(1 . the clumping factor is less than 10. For 
an expansion speed of 50 km/s and a terminal speed of 
2300 km/s we have e c i « 20%, i.e. for the model parame- 
ters used in this paper, the error on the clumping factor 
caused by applying the periodic boundary condition is less 
than 20%. The assumption of a single shell is not a limita- 
tion, because the above analysis can be repeated for any 
number of shells to obtain the same result. 



3.5. Implementation 

We have implemented the pseudo-planar periodic box 
technique in VH-1, a hydrodynamics programme (or hy- 
drocode) developed at the University of Virginia (Blondin, 
personal communication). This code solves the Eulerian 
equations of hydrodynamics by a Lagrangian remap (LR) 
method, which first solves the Lagrangian equations of hy- 
drodynamics and then maps the updated quantities back 
onto the Eulerian grid. A major challenge in Eulerian hy- 
drodynamics, even in their LR incarnation, is to obtain 
reliable estimates of time-averaged quantities such as pres- 
sure and density at zone interfaces. All quantities in VH- 
1 are zone-centred, i.e. VH-1 does not use a staggered 
mesh. To determine the time averages, a Riemann prob- 
lem (Zel'dovich & Raizer I1966|l is set up, and solved, at 
each zone interface. A parabolic interpolation is used to 
provide an accurate guess to the values on cither side of 
the interface. This combination of the parabolic interpola- 
tion and the Riemann solver is referred to as the piecewise 
parabolic method (PPM; Colella & Woodward I1984|> and 
is the heart of VH-1. The details of the implementation 
are described in Paper I, the obvious exception being that 
the radiative force need not be evaluated in the present 
work. 



4. Test problems 

4.1. Uniformly expanding sphere 

A straightforward problem to test the pseudo-planar ap- 
proach is that of an expanding sphere with spatially con- 
stant density and pressure, and a radially increasing ve- 
locity (see e.g. Blondin & Lufkin 1993). For an initial con- 
dition with density po, pressure po and velocity v = r/to 
this problem has the analytic solution 



P : 



P0[J 



P 



to 

Po\j 



37 



(17) 



This test problem was used by Blondin & Lufkin to il- 
lustrate the use of geometry corrections to minimise ad- 
vection errors near the origin of a curvilinear coordinate 
system. For the present simulation, we have used the ver- 
sion of VH-1 available from the North Carolina State 
University website. This version does not include the ge- 
ometry corrections. To avoid problems with the derivation 
of the internal energy from remapped quantities, Blondin 
& Lufkin removed the pressure gradient from the momen- 
tum equation. For the pseudo-planar method, this would 
also remove the geometrical source term from Eq. (0 , thus 
greatly reducing the stringency of our test. We therefore 
include the pressure gradient in our calculations. 

As a point of reference, we choose an initial condition 
that produces physical parameters typical for the outer 
wind of a hot star: 



tn 



10 5 s, po = 10- lb g/cm 



Po 



5xl0~ 4 dyn/cm 2 .(18) 



The spatial grid has 96 points and extends from 0.1 to 
3.1 i?*, where i?» = 19 Rq. We compare both a spherical 
model, calculated by solving Eqs. ©-Q, and a pseudo- 
planar model in a stationary box, calculated by solv- 
ing Eqs. ©-(|SJ|, with the exact analytical solution. The 
Courant number is 0.25. The boundary conditions were 
chosen so as to impose a zero gradient on density and pres- 
sure, and a constant gradient on the velocity. In the case of 
the pseudo-planar method, the scaled variables (e.g. p) are 
implicitly scaled back to their physical counterpart (p) be- 
fore imposing the boundary condition. This is important 
when comparing the two methods, as a difference in the 
quality of the boundary conditions would swamp the rela- 
tively small errors intrinsic to the method. Figure [21 shows 
that the pseudo-planar model performs marginally better 
than the spherical model. This is gratifying, as the test is 
unfavourable to the pseudo-planar method, which has to 
produce a flat physical density and pressure by evolving a 
curved rescaled density and pressure. 

4.2. A moving shock tube in planar geometry 

The Sod shock tube (Sod I1978f) has become a standard 
test for all hydrodynamics codes. It is a special case of 
the Riemann problem, describing the evolution of an ar- 
bitrary discontinuity (Zel'dovich & Raizer 1966). A hypo- 
thetical membrane separates two regions of uniform den- 
sity, pressure and velocity, where at least one of these 
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Fig. 2. Relative difference between the analytic solution 
and the computer simulation for the uniformly expanding 
sphere test problem. Upper panel: density. Lower panel: 
velocity. 



quantities is different on either side of the membrane. At 
t = 0, the membrane is removed and the jump discontinu- 
ity breaks up into a shock wave, a contact discontinuity 
and a rarefaction wave. In planar geometry, this problem 
has a semi-analytical solution (Sod 119751 Laney^)98). 
The standard problem presented by Sod, with the gas at 
rest at both sides of the membrane, a density ratio of 8 
and a pressure ratio of 10, is quite unchallenging for any 
modern hydrocode. The physical problems to which these 
codes are applied generally present more extreme condi- 
tions. 

In radiatively driven stellar winds, density contrasts 
are often orders of magnitude larger than in the standard 
Sod problem, while the gas moves over the grid at super- 
sonic speeds. We therefore propose a shock tube with a 
density ratio of 800 and a pressure ratio of 1000, with all of 
the gas initially moving at the same supersonic speed. By 
Galilean invariance, the solution should not depend on this 
initial speed (except for an obvious displacement). On one 
side of the membrane, we have pi = 100, pi = 100, ui = 60, 
on the other p r = 0.125, p r — 0.1, u r = 60. Furthermore, 
to highlight possible differences between forward and re- 
verse waves, we set two such shock tubes back to back, so 
that we have shocks, contact discontinuities and rarefac- 
tion waves running in both directions. The initial situa- 
tion is a layer of dense, high-pressure gas separated by two 
membranes from the surrounding sparse, low-pressure gas, 
with everything moving at the same speed. The original 
thickness of the dense layer is 0.2. 

The initial speed is more than fifty times the adiabatic 
sound speed (the value of the adiabatic constant in this 
problem is 1.4). This is a more challenging test than the 
standard Sod problem, but by no means more extreme 
than the conditions in a stellar wind. In Fig. we com- 




5.2 5.4 5.6 



Fig. 3. A VH-1 solution of the density for the moving 
shock tube problem. The exact solution (solid line) is com- 
pared to a VH-1 simulation with resolution Ax = 0.005 
(crosses) and Ax = 0.0025 (dots). 

pare the analytical solution (solid line) with the numerical 
solution using VH-1 in planar geometry, with a step size of 
0.005 (crosses) and a Courant number of 0.25. The results 
are quite dismal, particularly for the inward rarefaction, 
which shows a lot of unphysical substructure. Note also 
that the inward running shock is not at the correct posi- 
tion. 

This substructure is present for a wide range of code 
parameters. Increasing the resolution by a factor of two 
(dots) gives an unsubstantial improvement, as do varia- 
tions in the restriction on the Courant number (we re- 
peated the calculation with Courant numbers between 0.1 
and 0.6 and found similar results). The artefacts occur re- 
gardless of whether the internal energy is remapped sep- 
arately or derived from the remapped total energy (see 
below) . They also occur for a wide range of shock flatten- 
ing parameters (see below) 

In hot-star wind simulations, the internal energy is 
an ill-conditioned part of the total energy (Feldmeier 
1995J. We found that deriving the internal energy from 
the remapped total energy can result in artificial spikes of 
low temperature. To alleviate this problem one can remap 
both the internal and total energy, and derive the internal 
energy from the total energy only when the Mach number 
of the flow is sufficiently low (Blondin, personal commu- 
nication). We found that, though this fix is important in 
our radiatively driven calculations, it does not help with 
the artefacts shown in Fig. [3] . Indeed, the models shown 
includes this fix, which at M = 60 means that the internal 
energy is never determined using the total energy. When 
the total energy is used, the results are worse. 

We find that shock flattening (Colella & Woodward 
11984(1 does little to improve the results. This technique 
flattens the interpolating parabola in the vicinity of a 
shock. In the calculations presented here we have used 
the flattening parameters a/ 1 ) = 0.75, lu^ — 5.0 and 
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Fig. 4. Two ZEUS solutions of the moving shock tube 
problem. The exact solution (solid line) is compared to 
a ZEUS calculation using the Van Leer scheme (crosses) 
and a calculation using PPA (dots). 



Fig. 5. Two VH-1 solutions of the moving shock tube 
problem. The exact solution (solid line) is compared to a 
a standard VH-1 calculation using PPM (crosses) and a 
calculation using linear interpolation (dots). 



e = 0.33. Even with extreme shock flattening (a/ 1 ' = 0, 
a/ 2 ) = 10.0 and e = 0) the artefacts do not disappear, 
and the amount of smearing at the shocks is unaccept- 
able. Flattening is usually only applied in the Lagrangian 
hydro step. Extending flattening to the remap step doesn't 
improve the results substantially. If the interpolating func- 
tion is totally flattened in very grid point (i.e. by setting 
the flattening coefficient /, = 1 for all i), VH-1 reproduces 
the Godunov scheme. In this case, all features in the test 
problem are completely wiped out, because the scheme is 
too dissipative. By trial and error we find that setting the 
flattening coefficient fi = 0.025 for all i gives acceptable 
results. 

As a comparison, we ran the same simulation using 
ZEUS, a widely used MHD code developed by Stone & 
Norman (1992J). The numerical algorithms used to solve 
the Eulerian equations in ZEUS are quite different from 
those adopted in VH-1. It does not use a Lagrangian 
remap, but directly solves the Eulerian equations. The so- 
lution is split into two parts: a source step and an advection 
step. Although a staggered mesh is used, interpolations 
are still needed to determine the time-averaged values of 
variables at zone interfaces. ZEUS provides a second-order 
method (Van Leer) and a third-order method (PPA, for 
piecewise parabolic advection) . Note that PPA is only half 
of PPM, as it incorporates the parabolic interpolation, but 
not the Riemann solver. We have applied both the Van 
Leer scheme and PPA to the moving shock tube prob- 
lem. Figure 0] shows that the Van Leer scheme (crosses) 
does much better on this problem than VH-1, while PPA 
(dots) produces the same kind of unphysical substructure. 
The resolution used in this calculation is Ax — 0.0025. At 
a resolution of 0.005 (not shown on figure), the Van Leer 
scheme doesn't manage to capture the shock and the con- 
tact discontinuity, but doesn't produce any spurious fea- 
tures either. Given the differences between the two codes, 



it is quite surprising that the artefacts in VH-1 and ZEUS 
PPA are so similar. It strongly suggests that they are the 
result of the parabolic interpolation scheme used in both 
codes. To confirm this, we artificially lowered the order of 
VH-1 by setting the quadratic term in the interpolating 
function to zero. The results for the linear interpolation 
scheme are much better than for the higher order PPM 
(Fig. EJ) - The difference is due to errors in the remap step. 
This can be seen from the fact that using linear interpo- 
lation in the remap step only gives essentially the same 
results. 

The artificial substructure in this test problem is sim- 
ilar to the rarefaction shocks found by Falle (2002) in a 
ZEUS test calculation of a different Riemann problem, us- 
ing the Van Leer scheme. We have found, as expected, that 
VH-1 performs very well on Falle 's test problem and that 
using PPA in ZEUS also produces artefacts. 

Figure [B] shows the density and the isothermal sound 
speed for a VH-1 simulation using the pseudo-planar 
method. The pseudo-planar method, by its very nature, 
can only solve problems in a spherical geometry. Planar 
geometry was mimicked by taking a very large initial po- 
sition of the box. The box velocity was taken equal to 
the initial speed of the gas. The pseudo-planar method 
(crosses) is compared to the exact solution and to the 
ZEUS Van Leer simulation described above. It is clear 
that the pseudo-planar method performs very well on this 
test problem. This shows that the main stumbling block 
for VH-1 and ZEUS PPA is the supersonic velocity with 
which the features move over the grid. Indeed, any box ve- 
locity reducing this velocity to less extremely supersonic 
values gives adequate results. With minor modifications, 
VH-1 can also be used as a pure Lagrangian code. In this 
mode, the results are of the same quality as those of the 
periodic box technique, which corroborates the above con- 
clusion. 
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Fig. 6. Density and sound speed for a periodic box solu- 
tion of the moving shock tube problem. The periodic box 
method (crosses) is compared to the exact solution (solid 
line) and the ZEUS Van Leer method (dots). 

In summary, we can say that the combination of 
parabolic interpolation and highly supersonic advection 
leads to unphysical substructure and incorrect shock 
speeds. The results can be improved by removing the ad- 
vection from the scheme (e.g. by using the periodic box 
technique) or by avoiding parabolic interpolation. The les- 
son to be drawn from these experiments is that for prob- 
lems that are rapidly advected over the computational 
grid, higher order schemes such as PPM do not necessarily 
give more accurate results than lower-order schemes. 

5. Application to stellar winds 

5.1. Initial condition 

To apply the moving periodic box simulation we select 
as an initial condition a region from the snapshot shown 
in Fig. ^ The region is chosen so that it contains a suf- 
ficiently large number of shells. This is necessary to be 
able to evolve a model over a long period of time, as the 
evolution of the structure is largely determined by the 
fact that shells in the box collide and form denser shells 
as the box moves outward, counteracting their pressure- 
expansion. If only a single shell remains, the calculation 



becomes meaningless as there is obviously no opportu- 
nity for further collisions. The boundaries of the box were 
set to generate a periodic velocity. The density and pres- 
sure were then slightly modified near the inner bound- 
ary of the box, to ensure periodicity in the corresponding 
scaled variables and thereby avoid introducing any addi- 
tional discontinuities. In practice, the initial condition is 
the region r = 46.1 — 94 i?* on the snapshot. The scaled 
density is made periodic by introducing a linear "correc- 
tion" between the inner boundary and some r\ , so that the 
modified scaled density is the same at the left and right 
boundaries and the correction vanishes at v\ . The pressure 
is made periodic in the same manner. Points beyond t\ are 
not modified. In the simulations shown r\ = 49.2 

5.2. Comparison with radiatively driven model 

5.2.1. Stationary periodic box 

To better evaluate possible differences between a moving 
periodic box model and the standard radiatively driven 
model, we first apply the pseudo-planar equations in a 
stationary periodic box. We can compare the stationary 
periodic box model with the driven model only in a lim- 
ited region of the r, t plane. This is illustrated in Fig. [7| 
The wide box represents the domain in space and time 
covered by a radiatively driven model. The tall box rep- 
resents the domain covered by a stationary box model. In 
principle, the two models can be compared in the square 
region where they overlap. However, the two models use 
different boundary conditions, so the comparison is only 
meaningful in the region that depends exclusively on the 
initial condition and not on the boundary conditions. This 
region is the filled triangle in Fig.0 The upper edge of this 
triangle is in fact a C+ characteristic (sec e.g. Zel'dovich & 
Raizer 1966 ), along which a signal from the inner bound- 
ary condition is carried outward. It need not be a straight 
line. In the context of the present paper it can be viewed 
as the outer edge of an expanding shell and the upshot 
of the above is that we can only compare the box model 
with the driven model for shells that have not crossed the 
boundary of the box. 

Figure [S] shows a snapshot at 100 ksec, for a station- 
ary periodic box model (solid line) and a radiatively driven 
model (dashed line). We show only the region in r where 
the comparison is meaningful. It is clear that the two mod- 
els agree very closely. The small differences are most ob- 
vious in the velocity plot (upper panel) and appear to be 
due to the residual level of radiative driving beyond 45 i?» . 

5.2.2. Moving periodic box 

Let us now compare the results for a moving periodic box 
model with a radiatively driven model. Figure shows the 
domains in the r,t plane covered by the moving box model 
and the radiatively driven model. The region (filled) where 
the two can be meaningfully compared is somewhat larger 
than for the stationary box, but still limited, this time by 
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Fig. 7. Regions in the space-time plane covered by a full 
radiatively driven model and a stationary periodic box 
model. The dark triangle indicates the area in space-time 
covered by both models, i.e. the area where results of the 
models can be compared 
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Fig. 8. Snapshot of velocity, density and temperature at 
100 ksec for a radiatively driven model (solid line) and a 
stationary periodic box model (dashed line) 

both C + and C_ characteristics. Again, the meaningful 
region can be viewed as those shells that have not crossed 
the boundary of the box. 

Figure [Till shows a snapshot at 100 ksec, for a moving 
periodic box model (solid line) and a radiatively driven 
model (dashed line). It is clear that the agreement is less 
than perfect. In particular, the temperature of the mov- 
ing box model shows features that are not present in the 
driven calculation. Both models have broad regions of rar- 



efied hot gas (such as the one extending from 79 to 84.5 
i?*). This is gas that has been previously heated and has 
remained hot due to the inefficiency of radiative cooling 
at low densities. (We recall that we solve the energy equa- 
tion including radiative cooling, using the Raymond et al. 
|1976j cooling curve) . In addition to these broad warm re- 
gions, the periodic box model has narrow regions of gas 
heated by nearby shocks. The most conspicuous pair of 
such narrow, hot regions is centred around 90 These 
narrow regions are not present in the case of a stationary 
periodic box. They appear due to the Galilean transfor- 
mation which is the only difference between the moving 
and stationary periodic box calculations. These narrow 
regions are perfectly physical: the left region is heated by 
the forward shock at ~ 89 i?*, the right region by the 
reverse shock at 90.5 i?*. Their narrowness is explained 
by the fact that the shocks have only been able to heat 
the gas during a 100 ksec time interval, and is consistent 
with the velocity at which the gas flows out of the shock. 
It is in fact their absence in the radiatively driven calcu- 
lation that is disturbing. As gas passes through a shock, 
part of its kinetic energy is converted into internal energy. 
As it moves away from the shock, the shock-heated gas 
cools by emitting radiation. A useful expression for the 
length of the cooling zone was given by Feldmeier (Trj)95). 
Using his Eq. (A9) we find that the cooling length for even 
the weakest outer-wind shocks is well-resolved by the nu- 
merical grid. The cooling length for the relatively strong 
(X = 3.7) reverse shock at 90.5 i?* in Fig. rjU| is a huge 
400 RJ 

The absence of hot gas behind the shocks of the ra- 
diatively driven and stationary periodic box model is a 
manifestation of the "collapse" of cooling zones. This is a 
problem that has plagued all hydrodynamical simulations 
of hot-star winds that include energy balance. It has been 
attributed to a global thermal instability leading to an os- 
cillation of the width of the cooling zone ( Feldmeier 1 1 995[) 
and to advective diffusion fCooper I1994|l . Advective dif- 
fusion is caused by the combination of advection with the 
typical shape of the cooling curve. Between 10 5 and 10 7 K, 
this curve falls off as T" 1 / 2 . When a steep temperature 
feature is advected over the grid, it is inevitably smeared 
by numerical diffusion. Due to the T~ 1//2 -dependence of 
the cooling, the cold side of the smeared-out feature cools 
more rapidly than the warm side. This effectively steepens 
the feature again, but also makes it narrower. In this way, 
the numerical scheme takes a bite out of the cooling zone 
at every time-step, until it is eaten away completely. Our 
results suggest that, at least in the outer wind, advective 
diffusion is the cause of the collapse of the cooling zones. 

6. Results for an example application 

Let us now apply the periodic box technique to the prob- 
lem it was designed for: the evolution of stellar wind struc- 
ture out to very large distances. In Fig. fl^we show the 
evolution of the density contrast (density divided by mean 
density) as a greyscale image. During the 7 Msec interval 
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Fig. 9. Regions in the space-time plane covered by a full 
radiatively driven model and a moving periodic box model. 
The dark trapezium indicates the area in space-time cov- 
ered by both models, i.e. the area where results of the 
models can be compared 
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Fig. 10. Snapshot of velocity, density and temperature 
at 100 ksec for a radiatively driven model (solid line) and 
a moving periodic box model (dashed line) 



shown, the box moves from 50 to 1300 This figure 
shows the considerable conceptual advantage of the peri- 
odic box method. By following a fraction of the wind as 
it moves away from the star, the relevant physical mecha- 
nisms (pressure expansion and shell collisions) are imme- 
diately apparent. Similar greyscale images can be made 
for the other variables. 

The evolution of structure can be usefully described 
by a number of statistical quantities, such as the clump- 
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Fig. 11. Greyscale image of the density contrast (density 
divided by mean density) as a function of space and time. 
The intensity scale has been truncated to highlight the 
kinematics of the shells. 



ing factor and the velocity dispersion. These descriptors 
have been defined in Paper I, and involve the temporal 
average of variables such as the density at a given posi- 
tion in the wind. This is somewhat cumbersome in a mov- 
ing box model, because a given position in the wind does 
not correspond to a unique position within the moving 
box, and calculating a time average requires some awk- 
ward bookkeeping. It is more convenient to replace the 
temporal averages by spatial averages using the following 
ergodic approximation 



A v / p{x,t')dt' 



t-T 



x+X 



p(x' , t)dx' , 



(19) 



which holds if the statistical properties remain stationary 
over a time At = L/v\,. We can then calculate the clump- 
ing factor and the velocity dispersion, which are shown 
in Fig. El The velocity dispersion declines gradually to 
reach values barely above the sound speed (« 13 km/s), 
although the strongest shock in the final state still has a 
jump velocity of 25 km/s fFig. 11311 . The clumping factor 
has the oscillating behaviour typical of the competition 
between pressure expansion and shell collisions. Between 
200 and 1300 i?*, the clumping factor ranges from 2.5 to 
6. 

These simulations show that, within the ID as- 
sumption, the wind remains structured over huge dis- 
tances. Observational diagnostics such as radio emission 
(both thermal and non-thermal) are inevitably influenced. 
Inferred mass-loss rates, that scale as \ffc\ would be over- 
estimated by a factor of around two. It appears plausible, 
again within the approximations, that shocks survive out 
to the large distances needed to explain the non-thermal 
radio emission. 
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Fig. 12. Clumping factor and velocity dispersion as a 
function of radius for the box model. 
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Fig. 13. Snapshot of density, velocity and temperature at 
the end of the box simulation (7 Msec). 



7. Summary 

We have presented an efficient technique to study the 
evolution of instability-generated structure far away from 
the star. Because it follows a small part of the structure, 
high spatial resolution can be maintained with a relatively 
small number of points. Furthermore, the number of depth 
points required does not depend on the time-span of the 
simulation. The small number of points makes for smaller 
files and facilitates the analysis of the results. As the box 
follows the flow, the gas is advected over the grid at much 
lower speeds than in a traditional model. This makes the 
condition on the Courant number much less restrictive, 
resulting in a substantial reduction in computation time. 
The moving box also has a conceptual advantage, in that 
it makes it easier to visualise the physical processes that 
are happening in the simulation. 



We have shown that the high Mach numbers with 
which the gas moves over the grid in traditional models 
are a cause for concern, as it can cause unphysical struc- 
ture to appear, and physical structure to disappear. In 
this sense, the box model is not only faster, but also more 
accurate. It is therefore a useful tool to check whether a 
solution depends on the Galilean frame in which it is ob- 
tained. Using the box model, we have shown that the wind 
remains clumped out to 1300 i?» and that weak shocks re- 
main present. 

A key limitation of the present method is its restriction 
to just one-dimension (ID), with focus on the extensive 
structure in radius, but not accounting at all for likelihood 
that in 2D or 3D models, Rayleigh- Taylor or thin-shell in- 
stabilities would break up the assumed azimuthal coher- 
ence (Vishniac I1994JI . Future work should thus focus on 
extending these period box techniques to 2D or 3D, in con- 
junction with recent efforts to develop multi-dimensional 
radiation-hydrodynamics simulations of the initial forma- 
tion of structure arising from the line-deshadowing insta- 
bility (Dessart k Owocki I2TOI Gomez & Williams USHjl ■ 
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